1 Data Formatting and Scrubbing

1.1 Packages

1.2 Data Directories

1.3 Setting Subject IDs and Parameters

1.4 Read in Block and EV betas

# Schaefer Parcels (400)
betas.block.sch<-read.csv(file = paste0(data.path.betas,"betas_block_schaefer_long.csv")) %>% 
  dplyr::select(-X) %>%
  mutate(Liquid=factor(block, levels = c("BASE","JUICE","NEUT","SALT"), 
                       labels = c("Baseline","Juice","Neutral","Saltwater")))
betas.ev.sch<-read.csv(file = paste0(data.path.betas,"betas_ev_schaefer_long.csv")) %>% 
  dplyr::select(-X) %>%
  separate(col = condition, into = c("Liquid","Reward"), sep = "_", remove = FALSE) %>%
  mutate(tentplot = tent+1, AllvsNoReward=ifelse(Liquid=="Noreward", yes = "NoReward", no = "Reward"))

# Subcortical Parcels (19)
betas.block.sub<-read.csv(file = paste0(data.path.betas,"betas_block_subcortical_long.csv")) %>% 
  dplyr::select(-X) %>%
  mutate(Liquid=factor(block, levels = c("BASE","JUICE","NEUT","SALT"), 
                       labels = c("Baseline","Juice","Neutral","Saltwater")))
betas.ev.sub<-read.csv(file = paste0(data.path.betas,"betas_ev_subcortical_long.csv")) %>% 
  dplyr::select(-X) %>%
  separate(col = condition, into = c("Liquid","Reward"), sep = "_", remove = FALSE) %>%
  mutate(tentplot = tent+1, AllvsNoReward=ifelse(Liquid=="Noreward", yes = "NoReward", no = "Reward"))

1.5 Read in the behavioral & self-report data, format/average data based on relevant conditions

#data.beta.block<-read.csv(paste0(data.path.betas,"betas_block.csv")) %>% select(-X)
#data.beta.EV<-read.csv(paste0(data.path.betas,"betas_EV.csv")) %>% select(-X)
data.behavior <- read.delim(data.path.behavior)
data.runkey <- read.csv(runkey.path)
data.selfreport<-read.csv(selfreport.path)

# BASELINE subject data (filter using data pipeline from dplyr)
base <- filter(data.behavior,runType=="baseline") %>%
  dplyr::select(subID,runID,runType,trial,RT,ACC,switch_trial,task,congruency) %>%
  mutate(ERR=ifelse(ACC==0,1,0))

base.means <- base %>% group_by(subID) %>% 
  summarise(n=n(), meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
            meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
baseRT.means <- base %>% group_by(subID) %>% filter(ACC==1) %>%
  summarise(n=n(), meanRT = mean(RT), seRT = sd(RT)/sqrt(n-1))

# INCENTIVE subject data (filter using data pipeline from dplyr)
incentive <- data.behavior %>% filter(runType=="incentive" & !is.na(RT)) %>%
  dplyr::select(subID,runID,runType,trial,RT,ACC,switch_trial,task,congruency,Feedback,
         moneyreward,money = moneysymbol) %>%
  inner_join(data.runkey, by=c("subID","runID")) %>%
  mutate(critRT = critRT, 
         subRewarded = (RT<critRT & ACC==1)*1,
         ERR=ifelse(ACC==0,1,0),
         moneyCode = factor(moneyreward, levels = c(1,2,4), labels = c(-1,0,1)),
         liqCode = factor(liquid,levels=c("saltwater","neutral","juice"),
                                       labels=c(-1,0,1)),
         rankCode = factor(rank,levels=c(1,2,3),labels=c(1,0,-1)),
         orderCode = factor(rank,levels=c(1,2,3),labels=c(-1,0,1)),
         liqCodeJvN = factor(liquid,levels=c("saltwater","neutral","juice"),
                                       labels=c(0,-1,1)),
         liqCodeSvN = factor(liquid,levels=c("saltwater","neutral","juice"),
                                       labels=c(-1,1,0)),
         subID=as.factor(subID)) %>%
  group_by(subID) %>%
  mutate(zRT = scale(RT,center = TRUE, scale = TRUE))
incentive$subID<-as.numeric(levels(incentive$subID))[incentive$subID]
incentive$moneyCode<-as.numeric(levels(incentive$moneyCode)[incentive$moneyCode])
incentive$liqCode<-as.numeric(levels(incentive$liqCode)[incentive$liqCode])
incentive$rankCode<-as.numeric(levels(incentive$rankCode)[incentive$rankCode])
incentive$orderCode<-as.numeric(levels(incentive$orderCode)[incentive$orderCode])
incentive$liqCodeJvN<-as.numeric(levels(incentive$liqCodeJvN)[incentive$liqCodeJvN])
incentive$liqCodeSvN<-as.numeric(levels(incentive$liqCodeSvN)[incentive$liqCodeSvN])

1.6 Summarising the incentive data based on relevant conditions

# summarise means of incentive data for each subject, across all conditions
incentive.means <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID) %>% 
  summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE), 
            seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
            meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
            meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means <- incentive %>% group_by(subID) %>% filter(!subID %in% sub.exclude,ACC==1) %>%
  summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))

# summarise by subject and congruency across all conditions
incentive.means.cong <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID,congruency) %>% 
  summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE), 
            seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
            meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
            meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means.cong <- incentive %>% group_by(subID,congruency) %>% filter(!subID %in% sub.exclude,ACC==1) %>%
  summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))

# summarise by subject and switch cost across all conditions
incentive.means.switch <- incentive %>% filter(!subID %in% sub.exclude, trial!=1) %>% group_by(subID,switch_trial) %>% 
  summarise(n=n(), meanRR = mean(subRewarded, na.rm=TRUE), 
            seRR = sd(subRewarded, na.rm=TRUE)/sqrt(n-1),
            meanACC = mean(ACC, na.rm=TRUE), seACC = sd(ACC, na.rm=TRUE)/sqrt(n-1),
            meanERR = mean(ERR, na.rm=TRUE), seERR = sd(ERR,na.rm=TRUE)/sqrt(n-1))
incentiveRT.means.switch <- incentive %>% group_by(subID,switch_trial) %>% filter(!subID %in% sub.exclude,ACC==1,trial!=1) %>%
  summarise(n=n(), meanRT = mean(RT, na.rm=TRUE), seRT = sd(RT)/sqrt(n-1))

# summarise the means of the incentive data for each subject, grouped by condition
incentive9RT.means <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID, liquid, money) %>%
  filter(!is.na(subRewarded), ACC==1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),
                                                  labels=c(1,0,-1))))) %>%
  ungroup(subID)
incentive9.means = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded)) %>% group_by(subID, liquid, money) %>% 
  summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
  ungroup(subID) %>%
  left_join(incentive9RT.means, by = c("subID","money","Liquid","Reward","MoneyLabel","moneyCode","liqCode"))

# summarise the means of the incentive data for each subject, grouped by congruency and condition
incentive9RT.means.cong <- incentive %>% filter(!subID %in% sub.exclude) %>% group_by(subID, congruency,liquid, money) %>%
  filter(!is.na(subRewarded), ACC==1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneycode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqcode = as.numeric(as.character(factor(Liquid,levels=c("juice","neutral","saltwater"),
                                                  labels=c(1,0,-1))))) %>%
  ungroup(subID)
incentive9.means.cong = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded)) %>% 
  group_by(subID, congruency, liquid, money) %>% 
  summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
  ungroup(subID) 

# summarise the means of the incentive data for each subject, grouped by switch cost and condition
incentive9RT.means.switch <- incentive %>% filter(!subID %in% sub.exclude) %>% 
  group_by(subID, switch_trial,liquid,money) %>%
  filter(!is.na(subRewarded), ACC==1, trial!=1) %>% summarise(n=n(), meanRT = mean(RT)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneycode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqcode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),
                                                  labels=c(1,0,-1)))),
         switchcode = as.numeric(as.character(factor(switch_trial, levels = c("switch","repeat"), labels = c(1,0))))) %>%
  ungroup(subID)
incentive9.means.switch = incentive %>% filter(!subID %in% sub.exclude, !is.na(subRewarded), trial!=1) %>% 
  group_by(subID, switch_trial, liquid, money) %>% 
  summarise(n = n(), meanRR = mean(subRewarded), meanACC = mean(ACC), meanERR = mean(ERR)) %>%
  ungroup(liquid) %>% dplyr::rename(Liquid=liquid) %>%
  mutate(Liquid = factor(Liquid, levels = c("juice","neutral","saltwater"), labels = c("Juice","Neutral","Saltwater")),
         Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))),
         MoneyLabel = factor(x = money, levels = c("$","$$","$$$$"), labels = c("Low", "Med", "High")),
         moneyCode = as.numeric(as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(-1,0,1)))),
         liqCode = as.numeric(as.character(factor(Liquid,levels=c("Juice","Neutral","Saltwater"),labels=c(1,0,-1))))) %>%
  ungroup(subID) 

# summarise the self report data
# motivation ratings
selfreport.motivation <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
  dplyr::select(subID, juice1_motivetrial,juice2_motivetrial,juice4_motivetrial,
         neut1_motivetrial,neut2_motivetrial,neut4_motivetrial,
         salt1_motivetrial,salt2_motivetrial,salt4_motivetrial) %>%
  gather(key = condition, value = motive_rating, juice1_motivetrial:salt4_motivetrial) %>%
  separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
  dplyr::select(-misc) %>%
  mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
         condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
         condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
         condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
         condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
         condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
         condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
         condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
         condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
  separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
  mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))
# liking ratings
selfreport.liking <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
  dplyr::select(subID, juice1_liketrial,juice2_liketrial,juice4_liketrial,
         neut1_liketrial,neut2_liketrial,neut4_liketrial,
         salt1_liketrial,salt2_liketrial,salt4_liketrial) %>%
  gather(key = condition, value = liking_rating, juice1_liketrial:salt4_liketrial) %>%
  separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
  dplyr::select(-misc) %>%
  mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
         condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
         condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
         condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
         condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
         condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
         condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
         condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
         condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
  separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
  mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))
# performance ratings
selfreport.perform <- data.selfreport %>% filter(!subID %in% sub.exclude) %>%
  dplyr::select(subID, juice1_performtrial,juice2_performtrial,juice4_performtrial,
         neut1_performtrial,neut2_performtrial,neut4_performtrial,
         salt1_performtrial,salt2_performtrial,salt4_performtrial) %>%
  gather(key = condition, value = perform_rating, juice1_performtrial:salt4_performtrial) %>%
  separate(col = condition, into = c("condition","misc"), sep="_", remove = TRUE) %>%
  dplyr::select(-misc) %>%
  mutate(condition = replace(condition, list = which(condition=="juice1"), values = "Juice_$"),
         condition = replace(condition, list = which(condition=="juice2"), values = "Juice_$$"),
         condition = replace(condition, list = which(condition=="juice4"), values = "Juice_$$$$"),
         condition = replace(condition, list = which(condition=="neut1"), values = "Neutral_$"),
         condition = replace(condition, list = which(condition=="neut2"), values = "Neutral_$$"),
         condition = replace(condition, list = which(condition=="neut4"), values = "Neutral_$$$$"),
         condition = replace(condition, list = which(condition=="salt1"), values = "Saltwater_$"),
         condition = replace(condition, list = which(condition=="salt2"), values = "Saltwater_$$"),
         condition = replace(condition, list = which(condition=="salt4"), values = "Saltwater_$$$$")) %>%
  separate(col = condition, into = c("Liquid","money"), sep = "_") %>%
  mutate(Reward = as.character(factor(money,levels=c("$","$$","$$$$"),labels=c(1,2,4))))

1.7 Summarising Bundled Beta Values in dorsal Anterior Cingulate Cortex (Knot 2 & 3, 4-6 ms)

#'dACC' parcel ids from the incentive integration analyses
parcels.dACC<-c(107,110,108,311,312) 

# Looping over 'dACC' parcels and create data frame of combining effects from both parcels
for (pid in parcels.dACC) { # pid<-108
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    dplyr::summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.dACC[1])  {
    data.dACC<-tmp
  } else {
    data.dACC<- rbind(data.dACC,tmp)
  }
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values by Reward Rate across the parcels 107, 108, 110, 311, 312
data.dACC.sum <- data.dACC %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.block = mean(value.block),
            value.ev = mean(mean_value.ev),
            value.bundle.mean = mean(value.bundle), 
            meanRR = mean(meanRR),
            meanERR = mean(meanRR), 
            meanACC = mean(meanACC),
            meanRT = mean(meanRT), meanRTlog = log(meanRT),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")),
         mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code) %>%
  unite(col = "Incentive9", Liquid:Reward, sep="_", remove = FALSE)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.dACC.sum$motive_rating.normed[is.nan(data.dACC.sum$motive_rating.normed)]<-0

# Summarise by subject
data.dACC.sum.sub <- data.dACC %>% group_by(subID) %>%
  summarise(value.block = mean(value.block),
            value.ev = mean(mean_value.ev),
            value.bundle.mean = mean(value.bundle), 
            meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

1.8 Summarising Bundled Beta Values in Ventromedial PFC

#'vmpfc' parcel ids from the incentive integration analyses
#parcels.vmpfc<-c(169,379) 
parcels.vmpfc<-c(168, 169, 173, 174, 379, 380, 381, 382)

# Looping over 'vmpfc' parcels and create data frame of combining effects from both parcels
for (pid in parcels.vmpfc) { # pid<-169
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.vmpfc[1])  {
    data.vmpfc<-tmp
  } else {
    data.vmpfc<- rbind(data.vmpfc,tmp)
  }
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values across the vmpfc parcels 169 and 379
data.vmpfc.sum <- data.vmpfc %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.vmpfc.sum$motive_rating.normed[is.nan(data.vmpfc.sum$motive_rating.normed)]<-0

# Summarise by subject
data.vmpfc.sum.sub <- data.vmpfc %>% group_by(subID) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

1.9 Summarising Bundled Beta Values in Nucleus Accumbens

# 'accumbens' parcel ids from the incentive integration analyses
parcels.accumbens<-c(350,351) 

# Looping over 'vmpfc' parcels and create data frame of combining effects from both parcels
for (pid in parcels.accumbens) { # pid<-350
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.accumbens[1])  {
    data.accumbens<-tmp
  } else {
    data.accumbens<- rbind(data.accumbens,tmp)
  }
  
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values across the accumbens parcels 350 and 351
data.accumbens.sum <- data.accumbens %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.accumbens.sum$motive_rating.normed[is.nan(data.accumbens.sum$motive_rating.normed)]<-0
  

# Summarise by subject
data.accumbens.sum.sub <- data.accumbens %>% group_by(subID) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

1.10 Summarising Bundled Beta Values in Caudate

# caudate parcel ids from the incentive integration analyses
parcels.caudate<-c(355,356) 

# Looping over 'caudate' parcels and create data frame of combining effects from both parcels
for (pid in parcels.caudate) { # pid<-355
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.caudate[1])  {
    data.caudate<-tmp
  } else {
    data.caudate<- rbind(data.caudate,tmp)
  }
  
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values across the accumbens parcels 350 and 351
data.caudate.sum <- data.caudate %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
mutate(mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.caudate.sum$motive_rating.normed[is.nan(data.caudate.sum$motive_rating.normed)]<-0


# Summarise by subject
data.caudate.sum.sub <- data.caudate %>% group_by(subID) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

1.11 Summarising Bundled Beta Values in Putamen

# 'putamen' parcel ids from the incentive integration analyses
parcels.putamen<-c(361,362) # putamen

# Looping over 'putamen' parcels and create data frame of combining effects from both parcels
for (pid in parcels.putamen) { # pid<-361
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sub %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sub %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.putamen[1])  {
    data.putamen<-tmp
  } else {
    data.putamen<- rbind(data.putamen,tmp)
  }
  
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values across the accumbens parcels 350 and 351
data.putamen.sum <- data.putamen %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.putamen.sum$motive_rating.normed[is.nan(data.putamen.sum$motive_rating.normed)]<-0

# Summarise by subject
data.putamen.sum.sub <- data.putamen %>% group_by(subID) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

1.12 Summarising Bundled Beta Values in Anterior Insula

# 'anterior insula' parcel ids from the incentive integration analyses
parcels.insula<-c(99,101,143,306,340) # bilateral insula
#parcels.insula_left<-c(99,101,143) # left insula
#parcels.insula_right<-c(306,340) # right insula
#insula_hemispheres<-c("left","right")

# Looping over 'insula' parcels and create data frame of combining effects from both parcels
for (pid in parcels.insula) { # pid<-361
  # subset data for each parcel
  # blocked/sustained effects data 
  tmp.block <- betas.block.sch %>% filter(Liquid!="Baseline", parcel.id==pid) %>%
    dplyr::rename(value.block = value) %>%
    mutate(Liquid = as.character(Liquid))
  
  # bundled effects & behavioral data for cue related activity
  tmp <- betas.ev.sch %>% filter(condition!="Noreward", parcel.id==pid, 
                                    (tentplot>=TR_cue[1] & tentplot<=TR_cue[2])) %>%
    group_by(subID,Liquid,Reward) %>%
    summarise(mean_value = mean(value)) %>% dplyr::rename(mean_value.ev = mean_value) %>%
    mutate(liqCode = as.numeric(as.character(factor(Liquid,levels=c("Saltwater","Neutral","Juice"),labels=c(-1,0,1)))),
           moneyCode = as.numeric(as.character(factor(Reward,levels=c("1","2","4"),labels=c(-1,0,1))))) %>%
    left_join(y = tmp.block, by = c("subID","Liquid")) %>%
    mutate(value.bundle = mean_value.ev + value.block) %>%
    inner_join(y = incentive9.means, by = c("subID","Liquid","Reward","liqCode","moneyCode")) %>%
    inner_join(y = selfreport.motivation, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.liking, by = c("subID","Liquid","Reward","money")) %>%
    inner_join(y = selfreport.perform, by = c("subID","Liquid","Reward","money"))
  #tmp$liqCode <- as.numeric(levels(tmp$liqCode)[tmp$liqCode])
  #tmp$moneyCode <- as.numeric(levels(tmp$moneyCode)[tmp$moneyCode])
  
  if (pid == parcels.insula[1])  {
    data.insula<-tmp
  } else {
    data.insula<- rbind(data.insula,tmp)
  }
  
}
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
# Summarise the bundled values across the accumbens parcels 350 and 351
data.insula.sum <- data.insula %>% group_by(subID, Liquid, Reward) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(mot9Code = liqCode + moneyCode, 
         meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         meanRR.centered = scale(meanRR, center = TRUE, scale = FALSE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         motive_rating.centered = scale(motive_rating, center = TRUE, scale = FALSE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.centered = scale(value.bundle.mean, center = TRUE, scale = TRUE),
         value.bundle.mean.mot9 = value.bundle.mean*mot9Code,
         meanRR.mot9 = meanRR*mot9Code,
         motive_rating.mot9 = motive_rating*mot9Code)
`summarise()` regrouping output by 'subID', 'Liquid' (override with `.groups` argument)
data.insula.sum$motive_rating.normed[is.nan(data.insula.sum$motive_rating.normed)]<-0

# Summarise by subject
data.insula.sum.sub <- data.insula %>% group_by(subID) %>%
  summarise(value.bundle.mean = mean(value.bundle), meanRR = mean(meanRR), 
            meanERR = mean(meanRR), meanACC = mean(meanACC),
            motive_rating = mean(motive_rating), liking_rating = mean(liking_rating),
            liqCode = mean(liqCode), moneyCode = mean(moneyCode)) %>%
  #ungroup(Liquid,Reward) %>% group_by(subID) %>%
  mutate(meanRR.normed = scale(meanRR, center = TRUE, scale = TRUE),
         motive_rating.normed = scale(motive_rating, center = TRUE, scale = TRUE),
         value.bundle.mean.normed = scale(value.bundle.mean, center = TRUE, scale = TRUE))
`summarise()` ungrouping output (override with `.groups` argument)

2 Motivational Incentive Integration Effects: Reward Rate and Self-Reported Motivation Ratings

binom.test(x = 117, n = 288, p = .05, alternative = "two.sided")

    Exact binomial test

data:  117 and 288
number of successes = 117, number of trials = 288, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.05
95 percent confidence interval:
 0.3490173 0.4654323
sample estimates:
probability of success 
               0.40625 

2.1 Extended Data Figure 1-2: RT between baseline and incentive block

cond.RT <- data.behavior %>% dplyr::select(subID,runID,trial,condition=runType,RT,ACC,switch_trial) %>%
  filter(ACC==1) %>% mutate(RT_log = log(RT))

RT.sum=summarySEwithin2(data=cond.RT, measurevar = "RT", withinvars = c("condition"), idvar = "subID")
Automatically converting the following non-factors to factors: condition

Attaching package: 'data.table'
The following objects are masked from 'package:dplyr':

    between, first, last
The following object is masked from 'package:purrr':

    transpose
pandoc.table(RT.sum)

--------------------------------------------------------------
 condition     N      RT     RTNormed    sd      se      ci   
----------- ------- ------- ---------- ------- ------- -------
 baseline    3467    857.8    858.7     374.9   6.367   12.48 

 incentive   11443   641.7    641.4      262    2.449    4.8  
--------------------------------------------------------------
p.RT.0<-ggplot(data = RT.sum, aes(x=condition, y = RT, fill=condition)) +
  geom_bar(position=position_dodge(width=0.8), color="black",
           stat="identity", width=0.8) + 
  geom_errorbar(position=position_dodge(width=0.8),
                aes(ymin=RT-ci, ymax=RT+ci), width=.2) +
  xlab("Block") + ylab("Response Time (ms)") +
  labs(fill = "Block") +
  coord_cartesian(ylim=c(500,1000)) +
  scale_fill_brewer(palette="Set1") +
  theme_classic() + 
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
    axis.title=element_text(size=12,face = "bold"),
    axis.text=element_text(size=14),
    legend.title=element_blank(),
    legend.position="none",
    legend.box.background = element_rect(colour = "black"),
    strip.text.x = element_text(size = 12))
p.RT.0

2.2 Extended Data Paired T-test: RT between baseline and incentive block

# summarizing to have one datapoint per subject per condition.
cond.RT.sum <- cond.RT %>%group_by(subID,condition) %>%
  dplyr::summarise(meanRT = mean(RT), meanRT_log = mean(RT_log))
`summarise()` regrouping output by 'subID' (override with `.groups` argument)
# baseline vs incentive block
t.test(x = subset(cond.RT.sum, condition=="baseline")$meanRT_log,
       y = subset(cond.RT.sum, condition=="incentive")$meanRT_log,
       paired = TRUE)

    Paired t-test

data:  subset(cond.RT.sum, condition == "baseline")$meanRT_log and subset(cond.RT.sum, condition == "incentive")$meanRT_log
t = 16.672, df = 46, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2573245 0.3279937
sample estimates:
mean of the differences 
              0.2926591 
# calculate cohen's d
cohen.d(subset(cond.RT.sum, condition=="baseline")$meanRT_log,
        subset(cond.RT.sum, condition=="incentive")$meanRT_log,
        na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)

Cohen's d

d estimate: 1.571705 (large)
95 percent confidence interval:
   lower    upper 
1.291783 1.851628 

2.3 Extended Data Figure 1-2: Accuracy between baseline and incentive block

# Compute averages for baseline and incentive blocks across all subjects
cond.ACC <- data.behavior %>% dplyr::select(subID,runID,trial,condition=runType,RT,ACC,switch_trial)

ACC.sum=summarySEwithin2(data=cond.ACC, measurevar = "ACC", withinvars = c("condition"), idvar = "subID")
Automatically converting the following non-factors to factors: condition
p.ACC.0<-ggplot(data = ACC.sum, aes(x=condition, y = ACC, fill=condition)) +
  geom_bar(position=position_dodge(width=0.8), color="black",
           stat="identity", width = 0.8) + 
  geom_errorbar(position=position_dodge(width=0.8),
                aes(ymin=ACC-ci, ymax=ACC+ci), width=.2) +
  xlab("Block") + ylab("Accuracy") +
  labs(fill = "Block") +
  coord_cartesian(ylim=c(.7,1)) +
  scale_fill_brewer(palette="Set1") +
  theme_classic() + 
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
    axis.title=element_text(size=12,face = "bold"),
    axis.text=element_text(size=14),
    legend.title=element_blank(),
    legend.position="none",
    legend.box.background = element_rect(colour = "black"),
    strip.text.x = element_text(size = 12))
p.ACC.0

2.4 Extended Data Paired T-test: Accuracy between baseline and incentive block

# summarizing to have one datapoint per subject per condition.
cond.ACC.sum <- cond.ACC %>% group_by(subID,condition) %>%
  dplyr::summarise(meanACC = mean(ACC, na.rm=TRUE))
`summarise()` regrouping output by 'subID' (override with `.groups` argument)
# baseline vs incentive block
t.test(x = subset(cond.ACC.sum, condition=="baseline")$meanACC,
       y = subset(cond.ACC.sum, condition=="incentive")$meanACC,
       paired = TRUE)

    Paired t-test

data:  subset(cond.ACC.sum, condition == "baseline")$meanACC and subset(cond.ACC.sum, condition == "incentive")$meanACC
t = 7.81, df = 46, p-value = 5.645e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.05551917 0.09407461
sample estimates:
mean of the differences 
             0.07479689 
# calculate cohen's d
cohen.d(subset(cond.ACC.sum, condition=="baseline")$meanACC,
        subset(cond.ACC.sum, condition=="incentive")$meanACC,
        na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)

Cohen's d

d estimate: 1.111597 (large)
95 percent confidence interval:
    lower     upper 
0.7520448 1.4711495 

2.5 Extended Data Figure 1-2: RT switch costs between baseline and incentive block

cond.RT.switch <- cond.RT %>%
  group_by(subID,condition,taskswitch=switch_trial) %>% filter(trial!=1, ACC==1) %>%
  summarise(meanRT = mean(RT)) %>% spread(key = taskswitch, value = meanRT) %>%
  rename(taskswitch="switch",taskrepeat="repeat") %>%
  mutate(switchcost = taskswitch-taskrepeat)
`summarise()` regrouping output by 'subID', 'condition' (override with `.groups` argument)
SW.sum.0<-cond.RT.switch %>% group_by(condition) %>%
  summarise(SW.mean = mean(switchcost), SW.sd = sd(switchcost))
`summarise()` ungrouping output (override with `.groups` argument)
SW.sum<-summarySEwithin2(data=cond.RT.switch, measurevar = "switchcost", 
                        withinvars = c("condition"), idvar = "subID")
Automatically converting the following non-factors to factors: condition
pandoc.table(SW.sum)

------------------------------------------------------------------------
 condition   N    switchcost   switchcostNormed    sd      se      ci   
----------- ---- ------------ ------------------ ------- ------- -------
 baseline    47     52.95           52.95         47.16   6.879   13.85 

 incentive   47     22.19           22.19         47.16   6.879   13.85 
------------------------------------------------------------------------
p.SW.1<-ggplot(data = SW.sum, aes(x=condition, y = switchcost, fill=condition)) +
  geom_bar(position=position_dodge(width=0.8), color="black",
           stat="identity", width=0.8) + 
  geom_errorbar(position=position_dodge(width=0.8),
                aes(ymin=switchcost-ci, ymax=switchcost+ci), width=.2) +
  xlab("Block") + ylab("RT switch cost (ms)") +
  labs(fill = "Block") + scale_fill_brewer(palette="Set1") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
    axis.title=element_text(size=12,face = "bold"),
    axis.text=element_text(size=14),
    legend.position="none",
    legend.box.background = element_rect(colour = "black"),
    strip.text.x = element_text(size = 12))
p.SW.1

2.6 Extended Data Paired T-test: RT switch costs between baseline and incentive block

# baseline vs incentive block
t.test(x = subset(cond.RT.switch, condition=="baseline")$switchcost,
       y = subset(cond.RT.switch, condition=="incentive")$switchcost,
       paired = TRUE)

    Paired t-test

data:  subset(cond.RT.switch, condition == "baseline")$switchcost and subset(cond.RT.switch, condition == "incentive")$switchcost
t = 3.1615, df = 46, p-value = 0.002777
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 11.17366 50.33584
sample estimates:
mean of the differences 
               30.75475 
# calculate cohen's d
cohen.d(subset(cond.RT.switch, condition=="baseline")$switchcost,
        subset(cond.RT.switch, condition=="incentive")$switchcost,
        na.rm = TRUE, pooled = TRUE, paired = TRUE, hedges.correction = FALSE)

Cohen's d

d estimate: 0.5636807 (medium)
95 percent confidence interval:
    lower     upper 
0.1824817 0.9448797 

2.7 Plot RT by switch and non-switch trials

RT.switch.sum=summarySEwithin2(data=incentive9RT.means.switch, measurevar = "meanRT",
                        withinvars = c("Liquid","money","switch_trial"), idvar = "subID")
Automatically converting the following non-factors to factors: money, switch_trial
p.RT.switch.bar<-ggplot(data = RT.switch.sum, aes(x = money, y = meanRT, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
  facet_grid(.~switch_trial) + 
  xlab("Money") + ylab("RT (ms)") +
  labs(fill = "Liquid") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(600,700)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RT.switch.bar

m.SW.2<-lme4::lmer(formula = meanRT ~ switchcode * liqcode * moneycode + (1|subID), 
            data = incentive9RT.means.switch)
tab_model(m.SW.2)
  meanRT
Predictors Estimates CI p
(Intercept) 629.48 592.72 – 666.24 <0.001
switchcode 22.56 13.67 – 31.45 <0.001
liqcode -5.52 -13.22 – 2.18 0.160
moneycode -14.31 -22.00 – -6.61 <0.001
switchcode * liqcode 7.07 -3.82 – 17.96 0.203
switchcode * moneycode -3.36 -14.25 – 7.52 0.545
liqcode * moneycode 0.54 -8.89 – 9.97 0.910
(switchcode * liqcode) *
moneycode
10.23 -3.10 – 23.57 0.133
Random Effects
σ2 4258.48
τ00 subID 15709.48
ICC 0.79
N subID 46
Observations 828
Marginal R2 / Conditional R2 0.017 / 0.790

2.8 Extended Data Save Figure 1-2: Combine RT, Accuracy, and RT switch cost plots

p.figure2<-cowplot::plot_grid(p.RT.0, p.ACC.0, p.SW.1, labels = c("a","b","c"), align = "h", ncol = 3)
ggsave(filename = "Supp_Figure2.eps", plot = p.figure2, 
       device="eps", dpi = 300,
       path = figure.path, width = 9, height = 3, scale = 1)

2.9 Reward Rate by Money and Liquid Incentives

m.RR<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)

m.RR.JvN<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum, subset = Liquid!="Saltwater")

m.RR.SvN<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum, subset = Liquid!="Juice")

tab_model(m.RR,m.RR.JvN,m.RR.SvN,
          title = "Reward Rate by Experimental Factors",
          dv.labels = c("Reward Rate","Juice v. Neutral","Saltwater v. Neutral"),
          pred.labels = c("Intercept","Money", "Liquid","Money:Liquid"),
          show.re.var = TRUE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)
Reward Rate by Experimental Factors
  Reward Rate Juice v. Neutral Saltwater v. Neutral
Predictors Estimates Statistic p Estimates Statistic p Estimates Statistic p
Intercept 0.69 43.02 <0.001 0.70 37.76 <0.001 0.70 40.39 <0.001
Money 0.03 4.41 <0.001 0.02 2.63 0.009 0.02 2.12 0.034
Liquid 0.01 2.29 0.022 0.00 0.08 0.936 0.03 2.11 0.034
Money:Liquid -0.00 -0.21 0.830 0.00 0.27 0.786 -0.01 -0.42 0.673
Random Effects
σ2 0.01 0.01 0.01
τ00 0.01 subID 0.01 subID 0.01 subID
N 46 subID 46 subID 46 subID
Observations 414 276 276
Marginal R2 / Conditional R2 0.028 / 0.533 0.020 / 0.659 0.030 / 0.495

2.10 RT and Accuracy by Money and Liquid Incentives

m.RT<-lme4::lmer(formula = meanRT ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)
m.ACC<-lme4::lmer(formula = meanACC ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)

2.11 Extended Data Table 1-2: RT and Accuracy by Money and Liquid Incentives

tab_model(m.RT, m.ACC, 
          dv.labels = c("Response Time","Accuracy"),
          pred.labels = c("Intercept","Money", "Liquid","Money:Liquid"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)
  Response Time Accuracy
Predictors Estimates Statistic p Estimates Statistic p
Intercept 641.67 34.34 <0.001 0.86 88.02 <0.001
Money -16.25 -4.93 <0.001 0.01 1.36 0.175
Liquid -1.83 -0.56 0.579 0.01 1.95 0.051
Money:Liquid 5.14 1.27 0.202 0.01 1.81 0.071
N 46 subID 46 subID
Observations 414 414
Marginal R2 / Conditional R2 0.010 / 0.842 0.012 / 0.450

2.12 Extended Table 1-3: Multiple Regression Between Reward Rate and Motivation Ratings

m.motive<-lme4::lmer(formula = motive_rating ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)

m.RR.motive<-lme4::lmer(formula = meanRR ~ moneyCode * liqCode + motive_rating + (1|subID), 
                      data = data.dACC.sum)

tab_model(m.motive,m.RR, m.RR.motive,
          title = "Hierarchical Regression with Motivation Ratings",
          dv.labels = c("Motivation Ratings","Reward Rate", "Reward Rate"),
          pred.labels = c("Intercept","Money", "Liquid","Money:Liquid","Motivation Ratings"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)
Hierarchical Regression with Motivation Ratings
  Motivation Ratings Reward Rate Reward Rate
Predictors Estimates Statistic p Estimates Statistic p Estimates Statistic p
Intercept 4.93 42.09 <0.001 0.69 43.02 <0.001 0.61 23.71 <0.001
Money 0.69 9.31 <0.001 0.03 4.41 <0.001 0.01 2.15 0.031
Liquid 0.75 10.04 <0.001 0.01 2.29 0.022 0.00 0.04 0.967
Money:Liquid 0.02 0.18 0.858 -0.00 -0.21 0.830 -0.00 -0.26 0.795
Motivation Ratings 0.02 4.46 <0.001
N 46 subID 46 subID 46 subID
Observations 414 414 414
Marginal R2 / Conditional R2 0.259 / 0.431 0.028 / 0.533 0.059 / 0.549

2.13 Hierarchical Regression with Reward Rate and Motivation Ratings

anova(m.RR,m.RR.motive)
refitting model(s) with ML (instead of REML)
Data: data.dACC.sum
Models:
m.RR: meanRR ~ moneyCode * liqCode + (1 | subID)
m.RR.motive: meanRR ~ moneyCode * liqCode + motive_rating + (1 | subID)
            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
m.RR           6 -612.69 -588.53 312.34  -624.69                         
m.RR.motive    7 -630.38 -602.20 322.19  -644.38 19.698  1  9.072e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.14 Reward Rate by Motivation Ratings By Subject (Scatterplot)

# Reward Rate by Self Report Motivation Ratings by Subject
# RColorBrewer::brewer.pal(8, "Set1") # to get the hex colors for set 1
cor.RR.motive<-cor.test(x = data.dACC.sum.sub$motive_rating.normed, y = data.dACC.sum.sub$meanRR.normed)
Rsquared=round(cor.RR.motive$estimate^2,3)
textlab1 <- paste("R^2 == ", Rsquared)
textlab2 <- paste("t(44)= ",round(cor.RR.motive$statistic,3),", p = ",round(cor.RR.motive$p.value,2), sep="")

p.linecolor<-brewer.pal(n=8, name = 'Set1')
p.RR.motive.bysub<-ggplot(data = data.dACC.sum.sub, aes(x = motive_rating.normed, y = meanRR.normed)) +
  geom_point(size = 2) +
  geom_smooth(method="lm", se = TRUE, col = p.linecolor[2]) +
  annotate("text", x = .15, y = -2.0, label = textlab1, color=p.linecolor[2], size = 4, parse=TRUE, hjust=0) +
  annotate("text", x = .15, y = -2.3, label = textlab2, color=p.linecolor[2], size = 4, parse=FALSE, hjust=0) +
  xlab("Self-Report Motivation Ratings") + ylab("Reward Rate") +
  #ggtitle("Average Values by Subject") +
  scale_color_brewer(palette="Set1") +
  theme_classic() + 
  theme(panel.grid.minor=element_blank(),
        plot.background = element_rect(fill = "transparent",colour = NA),
        rect = element_rect(fill = "transparent"),
        axis.title=element_text(size=12,face = "bold"),
        axis.text=element_text(size=16),
        legend.position="none",
        strip.text.x = element_text(size = 16))
p.RR.motive.bysub
`geom_smooth()` using formula 'y ~ x'

2.15 Reward Rate by Motivation Ratings by Monetary Reward and Liquid (Scatterplot)

#p.linecolor<-brewer.pal(n=9, name = "Spectral")
p.RR.motive<-ggplot(data = data.dACC.sum, 
                            aes(x = motive_rating.normed, y = meanRR.normed, color=Liquid)) +
  geom_point(size = 2) +
  stat_smooth(method="glm", se = FALSE) +
  facet_grid(.~MoneyLabel) +
  #annotate("text", x = .35, y = -2.0, label = textlab1, color=p.linecolor[4], size = 5, parse=FALSE, hjust=0) +
  #annotate("text", x = .35, y = -2.3, label = textlab2, color=p.linecolor[4], size = 5, parse=FALSE, hjust=0) +
  labs(x = "Self Report Motivation Ratings", y = "Reward Rate", color = "Liquid") +
  scale_color_brewer(palette="Set2") +
  #scale_color_manual(labels = c("$", "$$","$$$$")) +
  theme_bw() + 
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 14)) 
p.RR.motive
`geom_smooth()` using formula 'y ~ x'

2.16 Motivation Ratings by Monetary Reward and Liquid (Barplot)

RR.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "motive_rating",
                        withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.RR.motive.bar<-ggplot(data = RR.sum, aes(x = MoneyLabel, y = motive_rating, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=motive_rating-se, ymax=motive_rating+se), width=.2) +
  xlab("Money") + ylab("Self-Report Motivation Ratings") +
  labs(fill = "Liquid") +
  #geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(1,7)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.motive.bar

2.17 Reward Rate by Monetary Reward and Liquid (Barplot)

RR.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanRR",
                        withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.RR.mot9<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
  xlab("Money") + ylab("Reward Rate") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.5,.75)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="none",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.mot9

2.18 Manuscript Save Figure 1

p.figure1<-cowplot::plot_grid(p.RR.mot9, p.RR.motive.bar, p.RR.motive.bysub,
                              label_size = 16, ncol = 3, rel_widths = c(1,1.2,1))
`geom_smooth()` using formula 'y ~ x'
p.figure1

#labels = c("B","C")
ggsave(filename = "Figure1.eps", plot = p.figure1, 
       device="eps", dpi = 300,
       path = figure.path, width = 12, height = 3, scale = 1.2)
Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
not supported on this device: reported only once per page

2.19 Extended Data Figure 1-4: RT and Accuracy by Monetary Reward and Liquid (Barplot)

RT.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanRT",
                        withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.RT.mot9<-ggplot(RT.sum, aes(x=MoneyLabel, y=meanRT, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
  xlab("Money") + ylab("Response Time (ms)") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(600,700)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="none",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RT.mot9

ACC.sum=summarySEwithin2(data=data.dACC.sum, measurevar = "meanACC",
                        withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.ACC.mot9<-ggplot(ACC.sum, aes(x=MoneyLabel, y=meanACC, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanACC-se, ymax=meanACC+se), width=.2) +
  xlab("Money") + ylab("Accuracy") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.75,1)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.ACC.mot9

2.20 Extended Data Save Figure 1-4

p.figure4_sup<-cowplot::plot_grid(p.RT.mot9, p.ACC.mot9,
                              label_size = 16, ncol = 2, rel_widths = c(1,1.3))
p.figure4_sup

ggsave(filename = "Supplemental_Figure1-4.eps", plot = p.figure4_sup, 
       device="eps", dpi = 300,
       path = figure.path, width = 7, height = 3, scale = 1.2)

3 Dorsal Anterior Cingulate Cortex (dACC) Encodes Both Monetary and Liquid Incentives in Bundled Beta Estimates

3.1 Block Beta Values by Liquid in dACC

m.dACC.RR.block<-lme4::lmer(formula = value.block ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)

3.3 Bundled Beta Values by Money in dACC

m.dACC.RR.bundle<-lme4::lmer(formula = value.bundle.mean ~ moneyCode * liqCode + (1|subID), 
                      data = data.dACC.sum)

3.4 Table of the Block, Event-Related,and Bundled Models

tab_model(m.dACC.RR.block,m.dACC.RR.cue, m.dACC.RR.bundle,
          #title = "Dorsal ACC Betas",
          dv.labels = c("Block Betas", "Event-Related Betas", "Bundled Betas"),
          pred.labels = c("Intercept","Money","Liquid","Money:Liquid"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)
  Block Betas Event-Related Betas Bundled Betas
Predictors Estimates Statistic p Estimates Statistic p Estimates Statistic p
Intercept -0.02 -1.51 0.131 0.17 9.04 <0.001 0.14 6.57 <0.001
Money 0.00 0.00 1.000 0.04 5.55 <0.001 0.04 4.66 <0.001
Liquid 0.02 4.58 <0.001 -0.00 -0.31 0.759 0.02 2.98 0.003
Money:Liquid 0.00 0.00 1.000 0.00 0.31 0.756 0.00 0.26 0.794
N 46 subID 46 subID 46 subID
Observations 414 414 414
Marginal R2 / Conditional R2 0.022 / 0.575 0.032 / 0.570 0.032 / 0.565

3.5 Manuscript Figure 2: Plotting dACC Blocked (Liquid) Effects & Cue (Money) Effects

dACC.block.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.block", 
                 withinvars = "Liquid", idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.figure2.block<-ggplot(data = dACC.block.sum, aes(x = Liquid, y = value.block, fill = Liquid)) +
  geom_bar(stat = "identity", color="black") +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=value.block-se, ymax=value.block+se), width=.2) +
  xlab("Liquid") + ylab("Block Beta Estimates") +
  labs(fill = "Liquid") +
  #geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  #coord_cartesian(ylim=c(.35,.8)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic()+
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="none",
        legend.title=element_text(size=14, face = "bold"),
        strip.text.x = element_text(size = 12))
p.figure2.block

dACC.ev.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.ev", 
                 withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.figure2.ev<-ggplot(data = dACC.ev.sum, aes(x = MoneyLabel, y = value.ev, fill = Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=value.ev-se, ymax=value.ev+se), width=.2) +
  xlab("Money") + ylab("Cue Beta Estimates") +
  labs(fill = "Liquid") +
  #geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  #coord_cartesian(ylim=c(.35,.8)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.figure2.ev

3.6 Manuscript Figure 2: Plotting dACC Bundled Effects

dACC.bundle.sum=summarySEwithin2(data = data.dACC.sum, measurevar = "value.bundle.mean", 
                 withinvars = c("Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: Liquid
p.figure2.bundle<-ggplot(data = dACC.bundle.sum, aes(x = MoneyLabel, y = value.bundle.mean, fill = Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=value.bundle.mean-se, ymax=value.bundle.mean+se), width=.2) +
  xlab("Money") + ylab("Bundled Beta Estimates") +
  labs(fill = "Liquid") +
  #geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  #coord_cartesian(ylim=c(.35,.8)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme_classic() +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.figure2.bundle

3.7 Manuscript Figure 2: Draw Tentzero function of trialwise beta estimates

#'dACC' parcel ids from the incentive integration analyses
parcels.dACC<-c(107,110,108,311,312) 
parcels.vmpfc<-c(168, 169, 173, 174, 379, 380, 381, 382)
# colors.TENT<-c("firebrick1","firebrick3","firebrick4",
#                "seagreen1","seagreen3","seagreen4",
#                "steelblue1","steelblue3","steelblue4")
colors.TENT<-c("#62c1a3","#3A9276","#1D493B",
               "#f98c61","#EC4909","#9D3106",
               "#8B9FCB","#4C69A9","#2C3D63")

data.tentplot<-betas.ev.sch %>% filter(condition!="Noreward", parcel.id %in% parcels.dACC) %>%
  mutate(time=tentplot*2) %>%
  group_by(time,Liquid,Reward,condition) %>%
  #group_by(time,Reward) %>%
  dplyr::summarise(mean_value = mean(value), .groups = "keep") 

# add zero padding to beginning and end of tent functions

data.tentplot2<-data.tentplot %>% #dplyr::filter(time<=4) %>% 
  ungroup() %>%
  add_row(time = c(rep(0,9),rep(24,9)), Liquid = data.tentplot$Liquid[1:18],
          Reward = data.tentplot$Reward[1:18], condition = data.tentplot$condition[1:18],
          mean_value=0) %>%
  mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")))
  
# data.tentplot<-data.tentplot %>% #dplyr::filter(time<=4) %>% 
#   ungroup(time) %>%
#   mutate(time=c(rep(0,9),rep(24,9)), mean_value=0) %>%
#   bind_rows(data.tentplot) %>%
#   mutate(MoneyLabel = factor(x = Reward, levels = c(1,2,4), labels = c("Low", "Med", "High")))

p.trialbetas<-ggplot(data = data.tentplot2, mapping = aes(x = time, y = mean_value)) +
  geom_rect(mapping = aes(xmin=3.5,xmax=6.5,ymin=-.2,ymax=0.4), fill = "#e0e1dd", color = "#e0e1dd") +
  geom_line(aes(group=Reward, linetype = MoneyLabel, color = condition), size= .9) +
  #geom_point(aes(group=Reward, color = condition)) +
  facet_grid(.~Liquid) +
  scale_color_manual(values=colors.TENT, guide = FALSE) +
  #geom_text(x=5, y=.4, label="Cue", size = 6) +
  labs(x = "Time from Cue Onset (ms)", y = "Trial Beta Estimates", linetype = "Money")  +
  scale_x_continuous(breaks=c(seq(0,24,4)), labels=c(seq(0,24,4)), limits=c(0,24)) +
  scale_y_continuous(breaks=c(seq(-.2,.5,.2)), labels=c(seq(-.2,.5,.2)), limits=c(-.2,.5)) +
  theme_classic() +
  theme(panel.grid.minor=element_blank(),
        plot.background = element_rect(fill = "transparent",colour = NA),
        rect = element_rect(fill = "transparent"),
        #plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=15,face = "bold"),
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 16, face="bold"),
        legend.title=element_text(size=16,face="bold"),
        legend.text=element_text(size=12),
        legend.position="top")
p.trialbetas

3.8 Manuscript Save Figure 2

ggsave(filename = "Figure2_block.eps", plot = p.figure2.block, 
       device="eps", dpi = 300,
       path = figure.path, width = 5, height = 3, scale = 1.1)

ggsave(filename = "Figure2_ev.eps", plot = p.figure2.ev, 
       device="eps", dpi = 300,
       path = figure.path, width = 7, height = 4, scale = 1.1)

ggsave(filename = "Figure2_bundle.eps", plot = p.figure2.bundle, 
       device="eps", dpi = 300,
       path = figure.path, width = 7, height = 4, scale = 1.1)

ggsave(filename = "Figure2_trialbetas.eps", plot = p.trialbetas, 
       device="eps", dpi = 300,
       path = figure.path, width = 10, height = 4, scale = 1.1,
       bg = "transparent")

4 Reward Rate and Self-Report Motivation Ratings Predicted by Dorsal ACC Bundled Betas

4.1 Reward Rate predicted by Bundled Betas

m.dACC.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.dACC.sum)

m.dACC.RR.null<-update(m.dACC.RR, formula = ~ . -value.bundle.mean )
BF_dACC_RR = BIC(m.dACC.RR.null) - BIC(m.dACC.RR)

4.2 Self-Report Motivation Ratings Predicted by Dorsal ACC Bundled Betas

m.dACC.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.dACC.sum)

4.3 Plotting the Models

tab_model(m.dACC.RR, m.dACC.motive,
          title = "Reward Rate by Bundled dACC Betas",
          dv.labels = c("Reward Rate","Motivation Ratings"),
          pred.labels = c("Intercept","Incentive Conditions","dACC Bundled Betas"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)
Reward Rate by Bundled dACC Betas
  Reward Rate Motivation Ratings
Predictors Estimates CI Statistic p Estimates CI Statistic p
Intercept 0.68 0.65 – 0.71 40.93 <0.001 4.70 4.45 – 4.94 37.43 <0.001
Incentive Conditions 0.02 0.01 – 0.03 2.91 0.004 0.67 0.54 – 0.81 9.78 <0.001
dACC Bundled Betas 0.09 0.02 – 0.17 2.41 0.016 1.64 0.83 – 2.45 3.96 <0.001
N 46 subID 46 subID
Observations 414 414
Marginal R2 / Conditional R2 0.040 / 0.591 0.294 / 0.501

5 Dorsal ACC Effects on Reward Rate Mediated by Self-Report Motivation Ratings (Within-Subjects Mediation Analyses)

x = dACC /newline y = Reward Rate /newline m = Motivation Ratings newline

5.1 Format the Data for Within-Subjects Mediation Analysis

data.dACC.mediation <- data.dACC.sum %>% 
  dplyr::select(subID,meanRR,motive_rating,value.bundle.mean,mot9Code,liqCode,moneyCode)

5.2 Step 0: Motivation Conditions Predict Reward Rate (Y), dACC Betas, and Motivation Ratings

m.step0.RR<-lme4::lmer(formula = meanRR ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
m.step0.beta<-lme4::lmer(formula = value.bundle.mean ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
m.step0.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + (1+mot9Code|subID), data = data.dACC.sum)
tab_model(m.step0.RR, m.step0.beta, m.step0.motive,
          dv.labels = c("Reward Rate","dACC Betas", "Motivation Ratings"),
          pred.labels = c("Intercept","Incentive Conditions"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)
  Reward Rate dACC Betas Motivation Ratings
Predictors Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p
Intercept 0.69 0.66 – 0.73 43.02 <0.001 0.14 0.10 – 0.18 6.57 <0.001 4.93 4.70 – 5.16 42.09 <0.001
Incentive Conditions 0.02 0.01 – 0.03 3.36 0.001 0.03 0.01 – 0.04 4.04 <0.001 0.72 0.58 – 0.86 10.35 <0.001
N 46 subID 46 subID 46 subID
Observations 414 414 414
Marginal R2 / Conditional R2 0.025 / 0.597 0.031 / 0.613 0.259 / 0.492

5.3 Step 1: Reward Rate (Y) by dACC Betas (X)

m.dACC.step1<-lme4::lmer(formula = meanRR ~ value.bundle.mean + (1|subID), 
                          data = data.dACC.sum)
tab_model(m.dACC.step1)
  meanRR
Predictors Estimates CI p
(Intercept) 0.67 0.64 – 0.71 <0.001
value.bundle.mean 0.14 0.07 – 0.22 <0.001
Random Effects
σ2 0.01
τ00 subID 0.01
ICC 0.49
N subID 46
Observations 414
Marginal R2 / Conditional R2 0.035 / 0.510

5.4 Step 2: Motivation Ratings (M) Predicted by dACC Betas

m.dACC.step2<-lme4::lmer(formula = motive_rating ~ value.bundle.mean + (1|subID), 
                         data = data.dACC.sum)

5.5 Step 3: dACC Betas Predict Reward Rate Controling for Motivation Ratings

m.dACC.step3<-lme4::lmer(formula = meanRR ~ value.bundle.mean + motive_rating + (1|subID), 
                         data = data.dACC.sum)

5.6 Plot all 3 Models

tab_model(m.dACC.step1,m.dACC.step2,m.dACC.step3,
          title = "dACC Within-Subjects Mediation Models",
          dv.labels = c("Reward Rate","Motivation Ratings","Reward Rate"),
          #pred.labels = c("Intercept","Money", "Liquid","Motivation Ratings","dACC Bundled Betas"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE, show.ci = FALSE)
dACC Within-Subjects Mediation Models
  Reward Rate Motivation Ratings Reward Rate
Predictors Estimates Statistic p Estimates Statistic p Estimates Statistic p
(Intercept) 0.67 40.87 <0.001 4.53 34.68 <0.001 0.59 26.08 <0.001
value.bundle.mean 0.14 3.77 <0.001 2.77 5.98 <0.001 0.08 2.15 0.032
motive_rating 0.02 5.53 <0.001
N 46 subID 46 subID 46 subID
Observations 414 414 414
Marginal R2 / Conditional R2 0.035 / 0.510 0.102 / 0.231 0.071 / 0.545

5.7 Step 4: Run the mediation using bmlm package

data.dACC.mediation<-data.frame(data.dACC.mediation)
data.dACC.mediation<-isolate(d = data.dACC.mediation,by = "subID",
                             value = c("meanRR","motive_rating","value.bundle.mean","mot9Code"))

# Running the multilevel mediation on mean-centered data
set.seed(200)
dACC.mediation.results = mlm(d = data.dACC.mediation, id = "subID",
                       x = "value.bundle.mean_cw", m = "motive_rating_cw", y = "meanRR_cw",
                       iter = 10000)
Estimating model, please wait.

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000211 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 16.9229 seconds (Warm-up)
Chain 1:                13.053 seconds (Sampling)
Chain 1:                29.9759 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 8.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 15.6859 seconds (Warm-up)
Chain 2:                13.3309 seconds (Sampling)
Chain 2:                29.0168 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000106 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 15.7006 seconds (Warm-up)
Chain 3:                13.9962 seconds (Sampling)
Chain 3:                29.6968 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000113 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.13 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 17.3432 seconds (Warm-up)
Chain 4:                13.5487 seconds (Sampling)
Chain 4:                30.8919 seconds (Total)
Chain 4: 
# visualizing the path
mlm_path_plot(dACC.mediation.results, text = T)

# summarising within-subject effects
mlm.sum<-mlm_summary(dACC.mediation.results, level = .9)
mlm_pars_plot(dACC.mediation.results, type = "coef", pars = c( "cp", "c", "me","pme"))
No id variables; using all as measure variables
Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
Please use `group_by()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.
Warning: `summarise_()` is deprecated as of dplyr 0.7.0.
Please use `summarise()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

mlm_pars_plot(dACC.mediation.results, type = "hist", pars = c( "cp", "c", "me"))
No id variables; using all as measure variables
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# check for convergence of mcmc traces, requires bayesplot package
mcmc_trace(as.data.frame(dACC.mediation.results), pars = c("a","b","cp","c"))

mcmc_rank_hist(as.data.frame(dACC.mediation.results), pars = c("a","b","cp","c"))

# Create a pretty table of the parameters estimates
formattable(mlm.sum)
Parameter Mean SE Median 5% 95% n_eff Rhat
a 3.23 0.99 3.23 1.60 4.84 16780 1
b 0.02 0.00 0.02 0.01 0.03 16523 1
cp 0.07 0.05 0.07 -0.01 0.15 16338 1
me 0.05 0.03 0.05 0.01 0.10 8200 1
c 0.13 0.05 0.13 0.04 0.21 20909 1
pme 0.37 20.48 0.42 0.05 1.10 20428 1

5.8 Inverse mediation, with RR mediating effects of dACC on motivation ratings.

set.seed(500)
dACC.mediation.inverse.results = mlm(d = data.dACC.mediation, id = "subID",
                       x = "value.bundle.mean_cw", m = "meanRR_cw", y = "motive_rating_cw",
                       iter = 10000)
Estimating model, please wait.

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00011 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.1 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 12.9923 seconds (Warm-up)
Chain 1:                6.2606 seconds (Sampling)
Chain 1:                19.2529 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 8.8e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 11.9064 seconds (Warm-up)
Chain 2:                9.0943 seconds (Sampling)
Chain 2:                21.0007 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 9.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 12.1626 seconds (Warm-up)
Chain 3:                12.2249 seconds (Sampling)
Chain 3:                24.3875 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'bmlm' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000137 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.37 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 11.9763 seconds (Warm-up)
Chain 4:                10.1281 seconds (Sampling)
Chain 4:                22.1044 seconds (Total)
Chain 4: 
mlm_path_plot(dACC.mediation.inverse.results, text = T)

mlm.sum.inverse<-mlm_summary(dACC.mediation.inverse.results, level = .9)
mlm_pars_plot(dACC.mediation.inverse.results, type = "coef", pars = c( "cp", "c", "me","pme"))
No id variables; using all as measure variables

# check for convergence of mcmc traces, requires bayesplot package
mcmc_trace(as.data.frame(dACC.mediation.inverse.results), pars = c("a","b","cp","c"))

mcmc_rank_hist(as.data.frame(dACC.mediation.inverse.results), pars = c("a","b","cp","c"))

# Create pretty table of the parameter estimates
formattable(mlm.sum.inverse)
Parameter Mean SE Median 5% 95% n_eff Rhat
a 0.12 0.05 0.12 0.04 0.21 15979 1
b 3.45 0.77 3.45 2.20 4.71 16692 1
cp 2.82 0.92 2.82 1.31 4.34 8112 1
me 0.45 0.26 0.44 0.06 0.89 11173 1
c 3.27 0.95 3.27 1.71 4.85 8351 1
pme 0.14 0.24 0.13 0.02 0.30 17530 1

6 Between-Subjects correlations between Reward Rate, Motivation Ratings, and dACC

fit1<-lm(formula = meanRR.normed  ~ value.bundle.mean.normed, data = data.dACC.sum.sub)
fit2<-lm(formula =  motive_rating.normed ~ value.bundle.mean.normed, data = data.dACC.sum.sub)
fit3<-lm(formula =  meanRR.normed ~ value.bundle.mean.normed + motive_rating.normed, data = data.dACC.sum.sub)

tab_model(fit1,fit2,fit3, 
          title = "Between-Subjects Effects",
          dv.labels = c("Reward Rate","Motivation Ratings","Reward Rate"),
          pred.labels = c("Intercept","DACC Betas","Motivation Ratings"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)
Between-Subjects Effects
  Reward Rate Motivation Ratings Reward Rate
Predictors Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p
Intercept 0.00 -0.29 – 0.29 0.00 1.000 0.00 -0.28 – 0.28 0.00 1.000 0.00 -0.29 – 0.29 0.00 1.000
DACC Betas 0.28 -0.02 – 0.57 1.90 0.064 0.33 0.04 – 0.62 2.32 0.025 0.23 -0.08 – 0.54 1.47 0.149
Motivation Ratings 0.15 -0.16 – 0.46 0.99 0.328
Observations 46 46 46
R2 / R2 adjusted 0.076 / 0.055 0.109 / 0.089 0.097 / 0.055

6.1 Manuscript Figure 3b: dACC Activation by Subject (Scatterplot)

6.2 Manuscript Save Figure 3b

p.dACC.bysub<-cowplot::plot_grid(p.dACC.RR.bysub, p.dACC.motive.bysub,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1.15, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.dACC.bysub

# Save the plots
ggsave(filename = "Figure3_dACC_bysubject.png", plot = p.dACC.bysub, 
        path = figure.path, device = "png", width = 4, height = 2, scale = 2.5, dpi = 300)

7 Dorsal ACC Selectively Encodes Subjective Motivational Value and Modulates Motivated Task Performance

7.1 Create table of mixed models of reward rate for forest plot

m.dACC.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.dACC.sum)
m.vmpfc.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.vmpfc.sum)
m.accumbens.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.accumbens.sum)
m.caudate.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.caudate.sum)
m.putamen.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.putamen.sum)
m.insula.RR<-lme4::lmer(formula = meanRR ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.insula.sum)

# Putting Models in a table 
tab_model(m.dACC.RR,m.vmpfc.RR,m.accumbens.RR,m.caudate.RR,m.putamen.RR,m.insula.RR,
          title = "Reward Rate Predicted by Bundled Betas",
          dv.labels = c("dACC","vmPFC","accumbens","caudate","putamen","anterior insula"),
          pred.labels = c("Intercept","Incentive Conditions","Bundled Betas"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)
Reward Rate Predicted by Bundled Betas
  dACC vmPFC accumbens caudate putamen anterior insula
Predictors Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p
Intercept 0.68 0.65 – 0.71 40.93 <0.001 0.70 0.66 – 0.73 42.70 <0.001 0.69 0.66 – 0.73 42.99 <0.001 0.69 0.66 – 0.72 42.93 <0.001 0.69 0.66 – 0.72 43.19 <0.001 0.69 0.66 – 0.73 42.87 <0.001
Incentive Conditions 0.02 0.01 – 0.03 2.91 0.004 0.02 0.01 – 0.03 3.34 0.001 0.02 0.01 – 0.03 3.37 0.001 0.02 0.01 – 0.03 3.26 0.001 0.02 0.01 – 0.03 3.10 0.002 0.02 0.01 – 0.03 3.21 0.001
Bundled Betas 0.09 0.02 – 0.17 2.41 0.016 0.01 -0.03 – 0.06 0.55 0.585 -0.00 -0.02 – 0.01 -0.58 0.562 0.03 -0.02 – 0.08 1.06 0.288 0.06 -0.01 – 0.13 1.68 0.093 0.02 -0.06 – 0.10 0.51 0.607
N 46 subID 46 subID 46 subID 46 subID 46 subID 46 subID
Observations 414 414 414 414 414 414
Marginal R2 / Conditional R2 0.040 / 0.591 0.026 / 0.598 0.026 / 0.596 0.028 / 0.593 0.031 / 0.594 0.026 / 0.595

7.2 Create table of mixed models of motivation ratings for forest plot

m.dACC.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.dACC.sum)
m.vmpfc.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.vmpfc.sum)
m.accumbens.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.accumbens.sum)
m.caudate.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.caudate.sum)
m.putamen.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.putamen.sum)
m.insula.motive<-lme4::lmer(formula = motive_rating ~ mot9Code + value.bundle.mean + (1+mot9Code|subID), 
                      data = data.insula.sum)

# Putting models into a table
tab_model(m.dACC.motive,m.vmpfc.motive,m.accumbens.motive,m.caudate.motive,m.putamen.motive,m.insula.motive,
          title = "Motivation Ratings Predicted by Bundled Betas",
          dv.labels = c("dACC","vmPFC","Accumbens","Caudate","Putamen","Anterior Insula"),
          pred.labels = c("Intercept","Incentive Conditions","Bundled Betas"),
          show.re.var = FALSE, show.stat = TRUE, show.icc = FALSE)
Motivation Ratings Predicted by Bundled Betas
  dACC vmPFC Accumbens Caudate Putamen Anterior Insula
Predictors Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p Estimates CI Statistic p
Intercept 4.70 4.45 – 4.94 37.43 <0.001 4.90 4.67 – 5.13 41.04 <0.001 4.93 4.70 – 5.16 42.06 <0.001 4.87 4.64 – 5.11 40.65 <0.001 4.87 4.64 – 5.09 42.32 <0.001 4.87 4.64 – 5.10 41.89 <0.001
Incentive Conditions 0.67 0.54 – 0.81 9.78 <0.001 0.72 0.58 – 0.86 10.38 <0.001 0.72 0.58 – 0.86 10.34 <0.001 0.71 0.57 – 0.84 10.27 <0.001 0.69 0.55 – 0.83 9.80 <0.001 0.68 0.55 – 0.82 9.85 <0.001
Bundled Betas 1.64 0.83 – 2.45 3.96 <0.001 -0.32 -0.84 – 0.21 -1.19 0.233 -0.02 -0.19 – 0.15 -0.22 0.825 0.74 0.11 – 1.37 2.29 0.022 1.29 0.51 – 2.07 3.22 0.001 1.36 0.43 – 2.29 2.87 0.004
N 46 subID 46 subID 46 subID 46 subID 46 subID 46 subID
Observations 414 414 414 414 414 414
Marginal R2 / Conditional R2 0.294 / 0.501 0.261 / 0.494 0.258 / 0.491 0.267 / 0.499 0.279 / 0.502 0.275 / 0.497

7.3 Combine Coefficients into Table for Forest Plot

# Reward Rate Coefs
dACC.coef.RR<-coef(summary(m.dACC.RR))
vmpfc.coef.RR<-coef(summary(m.vmpfc.RR))
accumbens.coef.RR<-coef(summary(m.accumbens.RR))
caudate.coef.RR<-coef(summary(m.caudate.RR))
putamen.coef.RR<-coef(summary(m.putamen.RR))
insula.coef.RR<-coef(summary(m.insula.RR))

# Motivation Rating Coefs
dACC.coef.motive<-coef(summary(m.dACC.motive))
vmpfc.coef.motive<-coef(summary(m.vmpfc.motive))
accumbens.coef.motive<-coef(summary(m.accumbens.motive))
caudate.coef.motive<-coef(summary(m.caudate.motive))
putamen.coef.motive<-coef(summary(m.putamen.motive))
insula.coef.motive<-coef(summary(m.insula.motive))

data.SV.names<-rep(c("dACC","vmPFC","NAc","Caudate","Putamen","Insula"),2)
data.SV.betas<-c(dACC.coef.RR["value.bundle.mean","Estimate"],
                 vmpfc.coef.RR["value.bundle.mean","Estimate"],
                 accumbens.coef.RR["value.bundle.mean","Estimate"],
                 caudate.coef.RR["value.bundle.mean","Estimate"],
                 putamen.coef.RR["value.bundle.mean","Estimate"],
                 insula.coef.RR["value.bundle.mean","Estimate"],
                 dACC.coef.motive["value.bundle.mean","Estimate"],
                 vmpfc.coef.motive["value.bundle.mean","Estimate"],
                 accumbens.coef.motive["value.bundle.mean","Estimate"],
                 caudate.coef.motive["value.bundle.mean","Estimate"],
                 putamen.coef.motive["value.bundle.mean","Estimate"],
                 insula.coef.motive["value.bundle.mean","Estimate"])
data.SV.std<-c(dACC.coef.RR["value.bundle.mean","Std. Error"],
               vmpfc.coef.RR["value.bundle.mean","Std. Error"],
               accumbens.coef.RR["value.bundle.mean","Std. Error"],
               caudate.coef.RR["value.bundle.mean","Std. Error"],
               putamen.coef.RR["value.bundle.mean","Std. Error"],
               insula.coef.RR["value.bundle.mean","Std. Error"],
               dACC.coef.motive["value.bundle.mean","Std. Error"],
               vmpfc.coef.motive["value.bundle.mean","Std. Error"],
               accumbens.coef.motive["value.bundle.mean","Std. Error"],
               caudate.coef.motive["value.bundle.mean","Std. Error"],
               putamen.coef.motive["value.bundle.mean","Std. Error"],
               insula.coef.motive["value.bundle.mean","Std. Error"])
data.SV.ci<-data.SV.std*1.96
data.SV.model<-c(rep("Reward Rate",6),rep("Motivation Ratings",6))
data.SV.number<-c(rep(1:6,2))

# Combine into a data frame
data.SV<-data.frame(names=data.SV.names,betas=data.SV.betas,se=data.SV.std,ci=data.SV.ci,Outcome=data.SV.model,number=data.SV.number)
data.SV$names<-factor(data.SV$names, levels = c("vmPFC","NAc","Caudate","Putamen","Insula","dACC"))

# Plot the Bundled Beta Estimates for Each of the ROIs
#p.linecolor<-brewer.pal(n=8, name = 'Set1')
p.bundled<-ggplot(data = data.SV, aes(x = betas, y = names, shape=Outcome, color=Outcome)) +
  geom_point(stat="identity", size = 3) +
  geom_errorbarh(aes(xmin=betas-ci, xmax=betas+ci, height = .1)) +
  geom_vline(xintercept=0, linetype="dashed", color="grey0", size = .7) +
  xlab("Bundled Beta Effect") + ylab("Regions of Interest") +
  #xlim(-0.05, .15) +
  facet_grid(.~Outcome, scales = "free_x")+
  scale_color_brewer(palette="Set1") +
  theme_classic() + 
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14, face = "bold"),
        axis.text=element_text(size=14),
        legend.position="none",
        strip.text.x = element_text(size = 14, face = "bold"))
p.bundled

7.4 Manuscript Save Figure 4: Forest Plot

p.figure4<-cowplot::plot_grid(p.bundled, 
                              labels = c("4"), label_size = 16)
ggsave(filename = "Figure4.eps", plot = p.figure4, 
       device="eps", dpi = 300,
       path = figure.path, width = 7, height = 5, scale = 1.1)

8 Exploratory Visulization and Plots (Not included in publication)

8.1 ROI Analyses: vmPFC

p.vmpfc<-cowplot::plot_grid(p.vmpfc.RR.bysub, p.vmpfc.motive.bysub,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1.15, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.vmpfc

8.1.1 Within-Subjects Scatterplots vmPFC

p.vmpfc.diff<-cowplot::plot_grid(p.vmpfc.MotDiff, p.vmpfc.MotDiff2,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.vmpfc.diff

8.2 ROI Analyses: Accumbens

p.accumbens<-cowplot::plot_grid(p.accumbens.RR.bysub, p.accumbens.motive.bysub,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.accumbens

8.3 ROI Analyses: Caudate

p.caudate<-cowplot::plot_grid(p.caudate.RR.bysub, p.caudate.motive.bysub,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1.15, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.caudate

8.4 Within-Subjects Scatterplots Caudate

p.caudate.diff<-cowplot::plot_grid(p.caudate.MotDiff, p.caudate.MotDiff2,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.caudate.diff

9 ROI Analyses: Putamen

p.putamen<-cowplot::plot_grid(p.putamen.RR.bysub, p.putamen.motive.bysub,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.putamen

9.1 Within-Subjects Scatterplots Putamen

p.putamen.diff<-cowplot::plot_grid(p.putamen.MotDiff, p.putamen.MotDiff2,
                              labels = c(" "," "), 
                              rel_widths = c(1, 1),
                              rel_heights = c(1, 1),
                              align="v", nrow = 1)
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
p.putamen.diff

9.2 Exploratory Visualization: Correlation Plot of bundled betas (by averaged across all conditions within a subject)

data.corr <- data.dACC.sum.sub %>% dplyr::select(subID, dACC = value.bundle.mean) %>%
  left_join(data.vmpfc.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc = value.bundle.mean) %>%
  left_join(data.accumbens.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens = value.bundle.mean) %>%
  left_join(data.caudate.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate = value.bundle.mean) %>%
  left_join(data.putamen.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate, putamen = value.bundle.mean) %>%
  left_join(data.insula.sum.sub, by = "subID") %>% dplyr::select(subID, dACC, vmpfc, accumbens, caudate, putamen,insula = value.bundle.mean) %>%
  dplyr::select(-subID)

ggcorrplot(corr = cor(data.corr), hc.order=TRUE, type = "lower",
           colors = c("#6D9EC1", "white", "#E46726"), 
           lab = TRUE,
           title = "Correlation Matrix of Bundled Beta Estimates Across ROIs", 
           outline.color = "black", show.diag = TRUE, 
           p.mat = cor_pmat(data.corr))

10 Exploratory Visualization: Congruency

RR.sum=summarySEwithin2(data=incentive.means.cong, measurevar = "meanRR",
                        withinvars = c("congruency"), idvar = "subID")
Automatically converting the following non-factors to factors: congruency
p.RR.cong<-ggplot(RR.sum, aes(x=congruency, y=meanRR)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
  xlab("Congruency") + ylab("Reward Rate") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.5,.75)) +
  #coord_cartesian(ylim=c(600,650)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.cong

RR.sum=summarySEwithin2(data=incentive9.means.cong, measurevar = "meanRR",
                        withinvars = c("congruency","Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: congruency
p.RR.mot9.cong<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
  xlab("Money") + ylab("Reward Rate") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.5,.8)) +
  facet_grid(~congruency)+
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.mot9.cong

ggsave(filename = "congruency.eps", plot = p.RR.cong, 
       device="eps", dpi = 300,
       path = figure.path, width = 3, height = 3, scale = 1.1)

ggsave(filename = "congruency_mot9.eps", plot = p.RR.mot9.cong, 
       device="eps", dpi = 300,
       path = figure.path, width = 6, height = 3, scale = 1.1)

11 Exploratory Visualization: Reward Rate switch costs

RR.sum=summarySEwithin2(data=incentive.means.switch, measurevar = "meanRR",
                        withinvars = c("switch_trial"), idvar = "subID")
Automatically converting the following non-factors to factors: switch_trial
p.RR.switch<-ggplot(RR.sum, aes(x=switch_trial, y=meanRR)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
  #xlab("Congruency") + 
  ylab("Reward Rate") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.5,.75)) +
  #coord_cartesian(ylim=c(600,650)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.switch

RR.sum=summarySEwithin2(data=incentive9.means.switch, measurevar = "meanRR",
                        withinvars = c("switch_trial","Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: switch_trial
p.RR.mot9.switch<-ggplot(RR.sum, aes(x=MoneyLabel, y=meanRR, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRR-se, ymax=meanRR+se), width=.2) +
  xlab("Money") + ylab("Reward Rate") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  coord_cartesian(ylim=c(.5,.8)) +
  facet_grid(~switch_trial)+
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RR.mot9.switch

ggsave(filename = "switchcostRR.eps", plot = p.RR.switch, 
       device="eps", dpi = 300,
       path = figure.path, width = 3, height = 3, scale = 1.1)

ggsave(filename = "switchcostRR_mot9.eps", plot = p.RR.mot9.switch, 
       device="eps", dpi = 300,
       path = figure.path, width = 6, height = 3, scale = 1.1)

12 Exploratory Visualization: RT Switch Costs

RT.sum=summarySEwithin2(data=incentiveRT.means.switch, measurevar = "meanRT",
                        withinvars = c("switch_trial"), idvar = "subID")
Automatically converting the following non-factors to factors: switch_trial
p.RT.switch<-ggplot(RT.sum, aes(x=switch_trial, y=meanRT)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
  #xlab("Congruency") + 
  ylab("Response Time") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  #coord_cartesian(ylim=c(.5,.75)) +
  coord_cartesian(ylim=c(600,700)) +
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RT.switch

RT.sum=summarySEwithin2(data=incentive9RT.means.switch, measurevar = "meanRT",
                        withinvars = c("switch_trial","Liquid","MoneyLabel"), idvar = "subID")
Automatically converting the following non-factors to factors: switch_trial
p.RT.mot9.switch<-ggplot(RT.sum, aes(x=MoneyLabel, y=meanRT, fill=Liquid)) +
  geom_bar(position=position_dodge(width=0.8), color="black", 
           stat="identity", width=0.8) +
  geom_errorbar(position=position_dodge(width=0.8), 
                aes(ymin=meanRT-se, ymax=meanRT+se), width=.2) +
  xlab("Money") + ylab("Response Time") +
  labs(fill = "Liquid") +
  geom_hline(yintercept = .4, linetype="dashed") +
  scale_fill_brewer(palette="Set2") +
  #coord_cartesian(ylim=c(.5,.8)) +
  coord_cartesian(ylim=c(550,700)) +
  facet_grid(~switch_trial)+
  #scale_fill_discrete(name="Monetary Reward") +
  theme(#plot.title=element_text(size=22,face="bold", vjust=2),
        axis.title=element_text(size=14,face = "bold"),
        axis.text=element_text(size=14),
        legend.position="right",
        legend.title=element_text(size=14, face = "bold"),
        legend.text = element_text(size=14),
        strip.text.x = element_text(size = 12))
p.RT.mot9.switch

#m.switch<-lmer(meanRT ~ switch_trial * liqCode * moneyCode + (1|subID), data = incentive9RT.means.switch)
#summary(m.switch)

ggsave(filename = "switchcostRT.eps", plot = p.RT.switch, 
       device="eps", dpi = 300,
       path = figure.path, width = 3, height = 3, scale = 1.1)

ggsave(filename = "switchcostRT_mot9.eps", plot = p.RT.mot9.switch, 
       device="eps", dpi = 300,
       path = figure.path, width = 6, height = 3, scale = 1.1)